Investigation of the incremental response of soils using a discrete element model 
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The incremental stress-strain relation of dense packings of polygons is investigated here by using 
molecular dynamics simulations. The comparison of the simulation results to the continuous theories 
is performed using explicit expressions for the averaged stress and strain over a representative 
volume element. The discussion of the incremental response raises two important questions of soil 
deformation: Is the incrementally non-linear theory appropriate to describe the soil mechanical 
response? Does a purely elastic regime exists in the deformation of granular materials?. In both 
cases our answer will be " no" . The question of stability is also discussed in terms of the Hill condition 
of stability for non-associated materials. 
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I. INTRODUCTION 

For many years the study of the mechanical behavior of 
soils was developed in the framework of linear elasticity 
[1] and the Mohr-Coulomb failure criterion [2] However, 
since the start of the development of the non- linear con- 
stitutive relations in 1968 [3], a great variety of constitu- 
tive models describing different aspects of soils have been 
proposed [4]. A crucial question has been brought for- 
ward: What it the most appropriate constitutive model 
to interpret the experimental result, or to implement a 
finite element code? Or more precisely, why is the con- 
stitutive relation I am using better than that one of the 
fellow next lab? 

In the last years, the discrete element approach has 
been used as a tool to investigate the mechanical response 
of soils at the grain level [5] . Several average procedures 
have been proposed to define the stress [6-8] and the 
strain tensor [9,10] in terms of the contact forces and 
displacements at the individual grains. These methods 
have been used to perform a direct calculation of the in- 
cremental stress-strain relation of assemblies of disks [11] 
and spheres [12], without any a-priori hypothesis about 
the constitutive relation. Some of the results lead to 
the conclusion that the non-associated theory of elasto- 
plasticity is sufficient to describe the observed incremen- 
tal behavior [11]. However, some recent investigations 
using three-dimensional loading paths of complex loading 
histories seem to contradict these results [13,12]. Since 
the simple spherical geometries of the grains overestimate 
the role of rotations in realistic soils [13], it is interest- 
ing to evaluate the incremental response using arbitrarily 
shaped particles. 

In this paper we investigate the incremental response 
in the quasi-static deformation of dense assemblies of 
polygonal particles. The comparison of the numerical 
simulations with the constitutive theories is performed 
by introducing the concept of Representative Volume El- 
ement (RVE). This volume is chosen the smear out the 
strong fluctuations of the stress and the deformation in 
the granular assembly. In the averaging, each grain is 



regarded as a piece of continuum. By supposing that the 
stress and the strain of the grain are concentrated at the 
small regions of the contacts, we obtain expressions for 
the averaged stress and strain over the RVE, in terms of 
the contact forces, and the individual displacements and 
rotations of the grains. The details of this homogeniza- 
tion method are presented in Sec. II. A short review 
of the incremental, rate-independent stress-strain mod- 
els is presented in Sec. III. We make special emphasis 
in the classical Drucker-Prager elasto-plastic models and 
the recently elaborated theory of hypoplasticity. The 
details of the particle model are presented in Sec. IV. 
The interparticle forces include elasticity, viscous damp- 
ing and friction with the possibility of slip. The system 
is driven by applying stress controlled tests on a rectan- 
gular framework consisting of four walls. Some loading 
programs were implemented in Sec. V, in order to lead to 
four basic question on the incremental response of soils: 
1) The existence of tensorial zones in the incremental 
response, 2) the validity of the superposition principle 
and 3) the existence of a finite elastic regime and 4) the 
question of stability according to the Hill condition. 



II. HOMOGENIZATION 

The aim of this section is to calculate the macro- 
mechanical quantities, the stress and strain tensors, from 
micro-mechanical variables of the granular assembly such 
as contact forces, rotations and displacements of individ- 
ual grains. 

A particular feature of granular materials is that both 
the stress and the deformation gradient are very con- 
centrated in small regions around the contacts between 
the grains, so that they vary strongly on short dis- 
tances. The standard homogenization procedure smears 
out these fluctuations by averaging these quantities over 
a RVE. The diameter d of the RVE must be such that 
5 <C d <C D, where S is the characteristic diameter of 
the particles and D is the characteristic length of the 
continuous variables. 
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FIG. 1. Representative volume element (RVE). 



We use here this procedure to obtain the averages of 
the stress and the strain tensors over a RVE in granular 
materials, which will allow us to compare the molecular 
dynamics simulations to the constitutive theories. Wc 
regard stress and strain to be continuously distributed 
through the grains, but concentrated at the contacts. It 
is important to comment that this averaging procedure 
would not be appropriate to describe the structure of 
the chain forces or the shear band because typical vari- 
ations of the stress corresponds to few particle diame- 
ters. Different averaging procedures using coarse-grained 
functions [8], or cutting the space in slide-shaped areas 
[14,10], can deal with the question of how one can per- 
form averages, and at the same time maintain these fea- 
tures. 

We will calculate the averages around a point x n of the 
granular sample taking a RVE calculated as follows: at 
the initial configuration, we select the grains whose cen- 
ter of mass are less than d from xq. Then the RVE is 
taken as the volume V enclosed by the initial configura- 
tion of the grains. See Fig. 1. The diameter d is taken, 
so that the averaged quantities are not sensible to the 
increase of the diameter by one particle diameter. 



A. Micro-mechanical stress 

The Cauchy stress tensor is defined using the force act- 
ing on an area element situated on or within the grains. 
Let / be the force applied on a surface element a whose 
normal unit vector is n. Then the stress is defined as the 
tensor satisfying [1]: 



Cfej^fc = lim a ->ofj/a, 



(1) 



where the Einstein summation convention is used. In ab- 
sence of body forces, the equilibrium equations in every 
volume element lead to [1]: 



daij/dxi = 0. 



(2) 



We are going to calculate the average of the stress ten- 
sor a over the RVE: 
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V 



TdV. 



(3) 



Since there is no stress in the voids of the granular 
media, the averaged stress can be written as the sum of 
integrals taken over the particles 
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(JijdV, 



(4) 



where V a is the volume of the particle a and N is the 
number of particles of the RVE. Using the identity 



d(xitjkj) 



dx k dx k 
Eq. (2), and the Gauss theorem, Eq. (4) leads to 



(5) 
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d(xj<Jki, 
dx k 



dV a 



Xi<Jkjrikda. 

(6) 



The right hand side is the sum over the surface inte- 
grals of each grain. dV a represents the surface of the 
grain a and n is the unit normal vector to the surface 
element da. 

An important feature of granular materials is that the 
stress acting on each grain boundary is concentrated in 
the small regions near to the contact points. Then we can 
use the definition given in Eq. (1) to express this stress 
on particle a in terms of the contact forces by introducing 
Dirac delta functions: 



(3=1 



(7) 



where x af3 and f al3 are the position and the force at the 
contact f3, and N a is the number of contacts of the par- 
ticle a. Replacing Eq. (7) into Eq. (6), we obtain 
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<?iA = — 



V ^ 



X ,- 



(8) 



Now we decompose x 



■a p 



c a +£ a/3 where x a is the po- 



sition of the center of mass and £ a/3 is the branch vector, 
connecting the center of mass of the particle to the point 
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of application of the contact force. Imposing this decom- 
position in Eq. (8), and using the equilibrium equations 
at each particle / Q/3 = 0we have 



a/3 



(9) 



From the equilibrium equations of the torques 
is symmetric, i. e., 



S/jC^ '/f 3 — £j 13 f? 13 ) = one obtains that this tensor 



Oij — (Tji = 0. 



(10) 



Then, the eigenvalues of this matrix are always real. 
This property leads to some simplifications, as we will 
see later. 



B. Micro- mechanical strain 

In elasticity theory, the strain tensor is defined as the 
symmetric part of the average of the displacement gradi- 
ent with respect to the equilibrium configuration of the 
assembly. Using the law of conservation of energy, one 
can define the stress-strain relation in this theory [1]. In 
granular materials, it is not possible to define the strain in 
this sense, because any loading involves a certain amount 
of plastic deformation at the contacts, so that it is not 
possible to define the initial reference state to calculate 
the strain. Nevertheless, one can define a strain tensor 
in the incremental sense. This is defined as the average 
of the displacement tensor taken from the deformation 
during a certain time interval. 

At the micro-mechanical level, the deformation of the 
granular materials is given by a displacement field u = 
r' — r at each point of the assembly. Here r and r' are 
the positions of a material point before and after defor- 
mation. The average of the strain and rotational tensors 
are defined as: 



q = -(F-F t ). 



(11) 



(12) 



Here F T is the transpose of the deformation gradient F, 
which is defined as 



(13) 



Using the Gauss theorem, we transform it into an in- 
tegral over the surface of the RVE 



V 



UiUjda, 



dV 



(14) 



where dV is the boundary of the volume element. We 
express this as the sum over the boundary particles of 
the RVE 



, N b 



V 



-i JdV, 



(15) 



where Nb is the number of boundary particles. To go 
further it is convenient to make some approximations. 
First, the displacements of the grains during deformation 
can be considered rigid except for the small deformations 
near to the contact that can be neglected. Then, if the 
displacements are small in comparison to the size of the 
particles, we can write the displacement of the material 
points inside of particle a as: 



Ui(x) w Axf + e ijk A<f)j(x k - x k ), 



(16) 



where Ax a , A<f) a and x a are displacement, rotation and 
center of mass of the particle a which contains the ma- 
terial point x, and e^fe is the anti-symmetric unit tensor. 
Setting a parameterization for each surface of the bound- 
ary grains over the RVE, the deformation gradient can 
be explicitly calculated in terms of grain rotations and 
displacements by replacing Eq. (16) in Eq. (15). 

In the particular case of a bidimensional assembly of 
polygons, the boundary of the RVE is given by a graph 
{x~\..x~2, XN b +\ — x\} consisting of all the edges be- 
longing to the external contour of the RVE, as shown in 
Fig. 1. In this case, Eq. (15) can be transformed as a 
sum of integrals over the segments of this contour. 

1 N >> fXp+l 

Fij = ttJ2 ^ + e ^A^(x k - 4)]nfa, (17) 

where = e^k and the unit vector ft 13 is perpendicular 
to the segment x t3 x^ + K Here corresponds to the index 
of the boundary segment. Aaf 3 , Acf) 13 and x@ are displace- 
ment, rotation and center of mass of the particle which 
contains this segment. Finally, by using the parameteri- 
zation x — x 13 + s(x l3+1 — x 13 ), where (0 < s < 1), we can 
integrate Eq. (17) to obtain 



TP. . 
r%2 



(18) 



where N r ] = e j)k (x k +1 -x{ ) and t= (x^ 1 -x' 3 )/2-x a . 
The stress tensor can be calculated taking the symmetric 
part of this tensor using Eq. (11). Contrary to the strain 
tensor calculated for spherical particles [15], the individ- 
ual rotation of the particles appears in the calculation 
of this tensor. This is given by the fact that for non- 
spherical particles the branch vector is not perpendicular 
to the contact normal vector, so that eiki^Nj ^ 0. 
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III. INCREMENTAL THEORY 

Since the stress and the strain tensor are symmetric, 
it is advantageous to simplify the notation by defining 
these quantities as six-dimensional vectors: 



on 
CT22 
CT33 

V2CT31 
V2a 13 



and e 



en 

£22 
£33 

V2e 23 
V2ei3 



(19) 



The coefficient V2 allows us to preserve the metric in 
this transformation: a^a^ — <Jij<Jij ■ The relation be- 
tween these two vectors will be established in the general 
context of the rate-independent incremental constitutive 
relations. We will focus on two particular theoretical de- 
velopments: the theory of hypoplasticity and the elasto- 
plastic models. The similarities and differences of both 
formulations are presented in the framework of the incre- 
mental theory as follows. 



A. General framework 

In principle, the mechanical response of granular mate- 
rials can be described by a functional dependence of the 
stress a(t) at time t on the strain history {e(i')}o<t'<t- 
However, the mathematical description of this depen- 
dence turns out to be very complicated due to the 
non-linearity and irreversible behavior of these mate- 
rials. An incremental relation, relating the incremen- 
tal stress da(t) = a 1 (t)dt with the incremental strain 
de{t) = e'(t)dt and some internal variables \ account- 
ing for the deformation history, enables us to avoid these 
mathematical difficulties [16]. This incremental scheme 
is also useful to solve geotechnical problems, since the fi- 
nite element codes require that the constitutive relation 
be expressed incrementally. 

Due to the large number of existing incremental rela- 
tions, the necessity of a unified theoretical framework has 
been pointed out as an essential necessity to classify all 
the existing models [17] In general, the incremental stress 
is related to the incremental strain by the following func- 
tion: 



T x {de, da, dt). 



(20) 



Let's look at the special case where there is no rate 
dependence in the constitutive relation. This means that 
this relation is not influenced by the time required dur- 
ing any loading tests, as corresponds to the quasi-static 
approximation. In this case T is invariant with respect 
to dt, and Eq. (20) can be reduced to: 



In particular, the rate-independent condition implies 
that multiplying the loading time by a scalar A does not 
affect the incremental stress-strain relation: 



VA, Q x (\da) = \Q x (da) 



(22) 



This equation means that Q x is an homogeneous func- 
tion of degree one. In this case, the application of the 
Euler identity shows that Eq. (21) leads to 



de = M x (a)da, 



(23) 



where M x = dQ x / d{da) and a is the unitary vector defin- 
ing the direction of the incremental stress 



a = 



da 

W\ 



(24) 



de = Q x {da) 



(21) 



Eq. (23) represents the general expression for the rate- 
independent constitutive relation. In order to determine 
the dependence of M on a, one can either perform ex- 
periments by taking different loading directions, or pos- 
tulate explicit expressions based on a theoretical frame- 
work. The first approach will be considered in the next 
section, and the discussion about some existing theoret- 
ical models will be presented as follows. 



B. Elasto-plastic models 

The classical theory of elasto-plasticity has been es- 
tablished by Drucker and Prager in the context if metal 
plasticity [18]. Some extensions have been developed to 
describe soils, clays, rocks, concrete, etc. [2,19]. Here 
we present a short review of these developments in soil 
mechanics. 

When a granular sample subjected to a confining pres- 
sure is loaded in the axial direction, it displays a typical 
stress-strain response as shown in the left part of Fig. 
2. A continuous decrease of the stiffness (i.e. the slope 
of the stress-strain curve) is observed during the load- 
ing. If the sample is unloaded, an abrupt increase in 
the stiffness is observed, leaving an irreversible deforma- 
tion. One observes that if the stress is changed around 
some region below a a, which is called the yield point, 
the deformation is almost linear and reversible. The first 
postulate of the elasto-plastic theory establishes a stress 
region immediately below the yield point where only elas- 
tic deformations are possible. 

Postulate 1: For each stage of loading there ex- 
ists a yield surface, which encloses a finite region 
in the stress space where only reversible deforma- 
tions are possible. 

The simple Mohr-Coulomb model assumes that the on- 
set of plastic deformation occurs at failure [2]. However, 
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it has been experimentally shown that plastic deforma- 
tion develops before failure [20]. In order to provide a 
consistent description of these experimental results with 
the elasto-plastic theory, it is necessary to suppose that 
the yield function changes with the deformation process 
[20]. This condition is schematically shown in Fig. 2. 
Let suppose that the sample is loaded until it reaches 
the stress a a and then it is slightly unloaded. If the 
sample is reloaded, it is able to recover the same stress- 
strain relation of the monotonic loading once it reaches 
the yield point a a again. If one increases the load to the 
stress <jb, a new elastic response can be observed after 
a loading reversal. In the elasto-plasticity context, this 
result is interpreted by supposing that the elastic regime 
is expanded to a new stress region below the yield point 

OB- 

Postulate 2: The yield function remains when 
the deformations take place inside of the elastic 
regime, and it changes as the plastic deformation 
evolves. 

The transition from the clastic to the elasto-plastic 
response is extrapolated for more general deformations. 
Part (b) of Fig. 2 shows the evolution of the elastic re- 
gion when the yield point moves in the stress space from 
a a to (t~b- The essential goal in the elasto-plastic theory 
is to find the correct description of the evolution of the 
elastic regime with the deformation, which is called the 
hardening law. 

We will finally introduce the third basic assumption 
from elasto-plasticity theory: 

Postulate 3: The strain can be separated in an 
elastic (recoverable) and a plastic (unrecoverable) 
component: 



de = de e + di p , 



(25) 



The incremental elastic strain is linked to the incre- 
mental stress by introducing an elastic tensor as 



da = D{a)d~e e 



(26) 



To calculate the incremental plastic strain, we intro- 
duce the yield surface as 



A 



dg 
da. 



(29) 



where g is the so-called plastic potential function, fol- 
lowing the Drucker-Prager postulates it can be shown 
that g = f [18]. However, a considerable amount of ex- 
perimental evidence has shown that in soils the plastic 
deformation is not perpendicular to the yield surfaces 
[21]. It is necessary to introduce this plastic potential to 
extend the Drucker-Prager models to the so-called non- 
associated models. 

The parameter A of Eq. (29) can be obtained from 
the so-called consistence condition. This condition comes 
from the second postulate, which establishes that after 
the movement of the stress state from a a to o~b = a~A+da 
the elastic regime must be expanded so that df = 0, as 
shown in Part (b) of Fig. 2. Using the chain rule one 
obtains: 



At d f A , d f dK A V 

df= dc7 i dai + d- K d? j d ^ 



0. 



(30) 



Replacing Eq. (29) in Eq. (30), we obtain the param- 
eter A 



A = ( d/ dn dg t df ^ 
dnde p -daj d(Ji % ' 



(31) 



We define the vectors Nf — df/dai and N* = dgjdoi 
and the unit vectors cf> = N y /\Nv\ and ip = N f /\Nf \ as 
the flow direction and the yield direction. In addition, 
the plastic modulus is defined as 



h = 



1 



df dn dg 



|jV"»||JV/| dn de v dui 



(32) 



Replacing Eq. (31) in Eq. (29) and using Eq. (32) we 
obtain: 



de p = —<p ■ da tjj. 



(33) 



f(a, k) = 0, 



(27) 



where k is introduced as an internal variable, which de- 
scribes the evolution of the clastic regime with the defor- 
mation. From experimental evidence, it has been shown 
that this variable can be taken as a function of the cu- 
mulative plastic strain [2,19] 



e p = / \J de k de t 
Jo 



dt 



(28) 



When the stress state reaches the yield surface, the 
plastic deformation evolves. This is assumed to be de- 
rived from a scalar function of the stress as follows: 





Strain "] 

FIG. 2. Evolution of the elastic regime a) stress-strain re- 
lation b) elastic regime in the stress space. 
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Note that this equation has been calculated for the case 
that the stress increment goes outside of the yield sur- 
face. If the stress increment takes place inside the yield 
surface, the second Drucker-Prager postulate establishes 
that de p = 0. Thus, the generalization of Eq. (33) is 
given by 

dl p = ^(4>-da) (34) 

where (x) — x when x > and (x) = otherwise. Fi- 
nally, the total strain response can be obtained from Eqs. 
(25) and (34): 

de = D- 1 {a)da + ^-da) $ (35) 

From this equation one can distinguish two zones in the 
incremental stress space where the incremental relation 
is linear. They are the so-called tensorial zones defined 
by Darve [16]. The existence of two tensorial zones, and 
the continuous transition of the incremental response at 
their boundary, are essential features of the elasto-plastic 
models. 

Although the elasto-plastic theory has shown to work 
well for monotonically increasing loading, it has shown 
some deficiencies in the description of complex loading 
histories [22]. There is also an extensive body of exper- 
imental evidence that shows that the elastic regime is 
extremely small and that the transition from elastic to 
an elasto-plastic response is rather smooth [4]. 

The concept of bounding surface has been introduced 
to generalize the classical elasto-plastic concepts [23] . In 
this model, for any given state within the surface, a 
proper mapping rule associates a corresponding image 
stress point on this surface. A measure of the distance 
between the actual and the image stress points is used to 
specify the plastic modulus in terms of a plastic modulus 
at the image stress state. Besides the versatility of this 
theory to describe a wide range of materials, it has the 
advantage that the elastic regime can be considered as 
vanishingly small, and therefore used to give a realistic 
description of unbound granular soils. 

It is the author's opinion that the most striking as- 
pect of the bounding surface theory with vanishing elas- 
tic range is that it establishes a convergence point for two 
different approaches of constitutive modeling: the elasto- 
plastic approaches originated from the Drucker-Prager 
theory, and the more recently developed hypoplastic the- 
ories. 

C. Hypoplastic models 

In recent years, an alternative approach for modeling 
soil behavior has been proposed, which departs from the 
framework of the elasto-plastic theory [24-26]. The dis- 
tinctive features of this approach are: 



• The absence of the decomposition of strain 
in plastic and elastic components. 

• The statement of a non-linear dependence 
of the incremental response with the strain 
rate directions. 

The most general expression has been provided by the 
so-called second order incremental non-linear models [25] . 
A particular class of these models which has received spe- 
cial attention in recent times is provided by the theory 
of hypoplasticity [24,26]. A general outline of this the- 
ory was laid down by Kolymbas [24] , leading to different 
formulations, such as the K-hypoplasticity developed in 
Karlsruhe [27], and the CLoE-hypoplasticity originated 
in Grenoble [26]. In the hypoplasticity, the most general 
constitutive equation takes the following form: 

da = L(a, v)de + N(a, u)\de\ , (36) 

where L is a second order tensor and N is a vector, both 
depending on the current state of the material, the stress 
a and the void ratio v. Hypoplastic equations provide a 
simplified description of loose and dense unbound gran- 
ular materials. A reduced number of parameters are in- 
troduced, which are very easy to calibrate [28]. 

In the theory of hypoplasticity, the stress-strain rela- 
tion is established by means of an incremental non-linear 
relation without any recourse to yield or boundary sur- 
faces. This non-linearity is reflected in the fact that the 
relation between the incremental stress and the incre- 
mental strain given in Eq. (36) is always non-linear. The 
incremental non-linearity of the granular materials is still 
under discussion. Indeed, an important feature of the 
incremental non-linear constitutive models is that they 
break away from the superposition principle: 

da{dl\ + de2) 7^ da(dei) + da(de2), (37) 

which is equivalent to: 

de{dd x + da 2 ) ^ de{dai) + de{da 2 ) (38) 




FIG. 3. Smooth and stair-like stress paths and correspond- 
ing strain responses, p and q represent the pressure and the 
deviatoric stress, e and 7 are the volumetric and deviatoric 
strain components. 
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An important consequence of this feature is addressed 
by Kolymbas [29] and Darve [25]. They consider two 
stress paths; the first one is smooth and the second 
one results from the superposition of small deviations 
as shown in Fig. 3. The superposition principle estab- 
lishes that the strain response of the stair-like path con- 
verges to the response of the proportional loading in the 
limit of arbitrarily small deviations. More precisely, the 
strain deviations Ae and the steps of the stress incre- 
ments Act satisfy limAo-^o Ae = 0. For the hypoplastic 
equations, and in general for the incremental non-linear 
models, this condition is never satisfied. For incremental 
relations with tensorial zones, this principle is satisfied 
whenever the increments take place inside one of these 
tensorial zones. It should be added that there is no ex- 
perimental evidence disproving or confirming this rather 
questionable superposition principle. 

IV. DISCRETE MODEL 

We present here a two-dimensional discrete element 
model which will be used to investigate the incremen- 
tal response of granular materials. This model consists 
of randomly generated convex polygons, which interact 
via contact forces. There are some limitations in the 
use of such a two-dimensional code to model physical 
phenomena that are three-dimensional in nature. These 
limitations have to be kept in mind in the interpretation 
of the results and its comparison with the experimental 
data. In order to give a three-dimensional picture of this 
model, one can consider the polygons as a collection of 
prismatic bodies with randomly-shaped polygonal basis. 
Alternatively, one could consider the polygons as three- 
dimensional grains whose centers of mass all move in the 
same plane. It is the author's opinion that it is more 
sensible to consider this model as an idealized granular 
material that can be used to check the constitutive mod- 
els. 

The details of the particle generation, the contact 
forces, the boundary conditions and the molecular dy- 
namics simulations are presented in this section. 

A. Generation of polygons 

The polygons representing the particles in this model 
are generated by using the method of Voronoi tessella- 
tion [30]. This methods is schematically shown in Fig. 4: 
First, a regular square lattice of side i is created. Then, 
we choose a random point in each cell of the rectangular 
grid. Then, each polygon is constructed assigning to each 
point that part of the plane that is nearer to it than to 
any other point. The details of the construction of the 
Voronoi cells can be found in the literature [31,32]. 




FIG. 4. Voronoi construction used to generate the convex 
polygons. The dots indicate the point used in the tessellation. 
Periodic boundary conditions were used. 

Using the Euler theorem, It has been shown analyt- 
ically that the mean number of edges of this Voronoi 
construction must be 6 [32] . The number of edges of the 
polygons is distributed between 4 and 8 for 98.7% of the 
polygons. It is also found that the orientational distri- 
bution of edges is isotropic, and the distribution of areas 
of polygons is symmetric around its mean value £ 2 . The 
probabilistic distribution of areas follows approximately 
a Gaussian distribution with variance of 0.36£ 2 . 

B. Contact forces 

In order to calculate the forces, we assume that all the 
polygons have the same thickness L. The force between 
two polygons is written as F = fL and the mass of the 
polygons is M = mL. In reality, when two elastic bodies 
come into contact, they have a slight deformation in the 
contact region. In the calculation of the contact force we 
suppose that the polygons are rigid, but we allow them 
to overlap. Then, we calculate the force from this virtual 
overlap. 

The first step for the calculation of the contact force is 
the definition of the line representing the flattened con- 
tact surface between the two bodies in contact. This is 
defined from the contact points resulting from the in- 
tersection of the edges of the overlapping polygons. In 
most cases, we have two contact points, as shown in the 
left of Fig. 5. In such a case, the contact line is de- 
fined by the vector C = C\C 2 connecting these two in- 
tersection points. In some pathological cases, the inter- 
section of the polygons leads to four or six contact points. 
In these cases, we define the contact line by the vector 
C = ChCf 2 + <Wl or C = C^f 2 + <Wt + CsCa, respec- 
tively. This choice guarantees a continuous change of the 
contact line, and therefore of the contact forces, during 
the evolution of the contact. 
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FIG. 5. Contact points d before (left) and after the for- 
mation of a pathological contact (right). The vector denotes 
the contact line, t represents the time step. 

The contact force is separated as f c = f e + f v , 
where f e and f v are the elastic and viscous contribu- 
tion. The elastic part of the contact force is decomposed 
as f e — .fn">i c + ftt°- The calculation of these compo- 
nents is explained below. The unit tangential vector is 
defined as i c — C/\C\, and the normal unit vector n c is 
taken perpendicular to C. The point of application of 
the contact force is taken as the center of mass of the 
overlapping polygons. 

As opposed to the Hertz theory for round contacts, 
there is no exact way to calculate the normal force be- 
tween interacting polygons. An alternative method has 
been proposed in order to model this force [33]. In 
this method, the normal elastic force is calculated as 
f% = —k n A/L c where k n is the normal stiffness, A is 
the overlapping area and L c is a characteristic length of 
the polygon pair. Our choice of L c is 1/2(1/ i?^ + 1/Rj) 
where Ri and Rj are the radii of the circles of the same 
area as the polygons. This normalization is necessary to 
be consistent in the units of force [30]. 

In order to model the quasi-static friction force, we cal- 
culate the elastic tangential force using an extension of 
the method proposed by Cundall-Strack [5]. An elastic 
force / t e = -fcfAsf proportional to the elastic displace- 
ment is included at each contact. k t is the tangential 
stiffness. The elastic displacement Axt is calculated as 
the time integral of the tangential velocity of the contact 
during the time where the elastic condition |/ t e | < fif^ is 
satisfied. The sliding condition is imposed, keeping this 
force constant when |/ t e | = fif^. The straightforward cal- 
culation of this elastic displacement is given by the time 
integral starting at the beginning of the contact: 

Ax$ = [ vWMtfZ - \ft\)dt', (39) 
Jo 

where is the Heaviside step function and v% denotes 
the tangential component of the relative velocity v c at 
the contact: 

v° = Vi — Vj — Ui x £i + ujj x £j . (40) 



Here Vi is the linear velocity and t<3j is the angular ve- 
locity of the particles in contact. The branch vector £j 
connects the center of mass of particle i with the point 
of application of the contact force. 

Damping forces are included in order to allow for rapid 
relaxation during the preparation of the sample, and to 
reduce the acoustic waves produced during the loading. 
These forces are calculated as 

fc = ~m( ln v c n h c + lt vci c ), (41) 

being m = (1/rrii + 1/rrij)^ 1 the effective mass of the 
polygons in contact. h c and t c are the normal and tan- 
gential unit vectors defined before, and 7„ and j t are the 
coefficients of viscosity. These forces introduce time de- 
pendent effects during the cyclic loading. However, we 
will show that these effects can be arbitrarily reduced by 
increasing the loading time, as corresponds to the quasi- 
static approximation. 

C. Molecular dynamics simulation 

The evolution of the position Xi and the orientation ipi 
of the i t h polygon is governed by the equations of motion: 

c b 

I i <Pi = 'Z$xf} + 1 E%xft. (42) 

c b 

Here mj and Ii are the mass and moment of inertia of 
the polygon i. The first summation goes over all particles 
in contact with this polygon. According to the previous 
section, these forces f c are given by 

f c = -(k n A/L c + j n mv c n )n c - {Ax c t + -f t mv L t )t c , 

(43) 

The second summation on the left hand of Eq. 43 
goes over all the vertices of the polygons in contact with 
the walls. This interaction is modeled by using a simple 
visco-elastic force. First, we allow the polygons to pene- 
trate the walls. Then, for each vertex of the polygon a 
inside of the walls we include a force 

/*• = -k n 5n - lb m a ^, (44) 

where S is the penetration length of the vertex, n is the 
unit normal vector to the wall, and v b is the relative ve- 
locity of the vertex with respect to the wall. 

We use a fifth-order Gear predictor-corrector method 
for solving the equation of motion [34]. This algorithm 
consists of three steps. The first step predicts position 
and velocity of the particles by means of a Taylor expan- 
sion. The second step calculates the forces as a function 
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of the predicted positions and velocities. The third step 
corrects the positions and velocities in order to optimize 
the stability of the algorithm. This method is much more 
efficient than the simple Euler approach or the Rungc- 
Kutta method, especially for problems where very high 
accuracy is a requirement. 

The parameters of the molecular dynamics simulations 
were adjusted according to the following criteria: 1) guar- 
antee the stability of the numerical solution, 2) optimize 
the time of the calculation, and 3) provide a reasonable 
agreement with the experimental data. 

There are many parameters in the molecular dynam- 
ics algorithm. Before choosing them, it is convenient to 
make a dimensional analysis. In this way, we can keep 
the scale invariance of the model and reduce the param- 
eters to a minimum of dimcnsionlcss constants. First, 
we introduce the following characteristic times of the 
simulations: the loading time t , the relaxation times 
t n = 1/7™, U = l/7t, tt — and the characteristic 
period of oscillation t s = \/k n /p£ 2 of the normal con- 
tact. 

Using the Buckingham Pi theorem [35], one can show 
that the strain response, or any other dimensionless vari- 
able measuring the response of the assembly during load- 
ing, depends only on the following dimensionless param- 
eters: a x = t n /t s , (X2 = t t /t s , a 3 = tb/t s , a-4 = to/t s , the 
ratio C = kt/k n between the stiffnesses, the friction coef- 
ficient /i and the ratio <Ji / k n between the stresses applied 
on the walls and the normal stiffness. 

The variables on will be called control parameters. 
They are chosen in order to satisfy the quasi-static ap- 
proximation, i.e. the response of the system does not 
depend on the loading time, but a change of these pa- 
rameters within this limit does not affect the strain re- 
sponse. a.2 — 0.1 and «2 = 0.5 were taken large enough 
to have a high dissipation, but not too large to keep the 
numerical stability of the method. 0:3 = 0.001 is chosen 
by optimizing the time of consolidation of the sample in 
the Subsec. IV D. The ratio a 4 = t /t s = 10000 was 
chosen large enough in order to avoid rate-dependence 
in the strain response, corresponding to the quasi-static 
approximation. Technically, this is performed by looking 
for the value of 0:4 such that a reduction of it by half 
makes a change of the stress-strain relation less than 5%. 

The parameters £ and fi can be considered as material 
parameters. They determine the constitutive response 
of the system, so they must be adjusted to the experi- 
mental data. In this study, we have adjusted them by 
comparing the simulation of biaxial tests of dense polyg- 
onal packings with the corresponding tests with dense 
Hostun sand [36]. First, the initial Young modulus of 
the material is linearly related to the value of normal 
stiffness of the contact. Thus, k n = l60MPa is chosen 
by fitting the initial slope of the stress-strain relation in 
the biaxial test. Then, the Poisson ratio depends on the 
ratio ( = k t /k n . Our choice k t — 52.8MPa gives an 



initial Poisson ratio of 0.07. This is obtained from the 
initial slope of the curve of volumetric strains versus ax- 
ial strain. The angles of friction and the dilatancy are 
increasing functions of the friction coefficient fi. Taking 
H = 0.25 yields angles of friction of 30 — 40 degrees and 
dilatancy angles of 10 — 20 degrees, which are similar to 
the experimental data of river sand [37] . 

D. Sample preparation 

The Voronoi construction presented above leads to 
samples with no porosity. This ideal case contrasts with 
realistic soils, where only porosities up to a certain value 
can be achieved. In this section, we present a method to 
create stable, irregular packings of polygons with a given 
porosity. 

The porosity can be defined using the concept of void 
ratio. This is defined as the ratio of the volume of the 
void space to the volume of the solid material. It can be 
calculated as: 



This is given in terms of the area enclosed by the floppy 
boundary V t , the sum of the areas of the polygons Vf and 
the sum of the overlapping areas between them V®. 

Of course, there is a maximal void ratio that can be 
achieved, because it is impossible to pack particles with 
an arbitrarily high porosity. The maximal void ratio v m 
can be detected by moving the walls until a certain void 
ratio is reached. We find a critical value, above which the 
particles can be arranged without touching. Since there 
is no contacts, the assembly cannot support a load, and 
even applying gravity will cause it to compactify. For a 
void ratio below this critical value, there will be parti- 
cle overlap, and the assembly is able to sustain a certain 
load. This maximal value is around 0.28. 

The void ratio can be decreased by reducing the vol- 
ume between the walls. The drawback of this method is 
that it leads to significant differences between the inner 
and outer parts of the boundary assembly, and it yields 
unrealistic overlaps between the particles, giving rise to 
enormous pressures. Alternatively, one can confine the 
polygons by applying gravity to the particles and on the 
walls. Particularly, homogeneous, isotropic assemblies, 
as shown in Fig. 6 can be generated by a gravitational 
field g — —gr where r is the vector connecting the center 
of mass of the assembly to the center of the polygon. 

When the sample is consolidated, repeated shear stress 
cycles are applied from the walls, superimposed to a con- 
fining pressure. The external load is imposed by applying 
a force [p c + q c sin(27rt/t )]W and [p c + q c cos(2Trt/to)]H 
on the horizontal and vertical walls, respectively. W and 
H are the width and the height of the sample. If we take 
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p c = 16kPa and q c < OAp c , we observe that the void 
ratio decreases as the number of cycles increases. Void 
ratios around 0.15 can be obtained. It is worth mention- 
ing that after some thousands of cycles the void ratio is 
still slowly decreasing, making it difficult to identify this 
minimal value. 

V. SIMULATION RESULTS 

In order to investigate different aspects of the incre- 
mental response some numerical simulations were per- 
formed. Different polygonal assemblies of 400 parti- 
cles were used in the calculations. The loading be- 
tween two stress states was controlled by applying forces 
[<t\ + {<t{ -<j{)r(t)]W and [a l 2 +(a^-(X l 2 )r(t)]H. A smooth 
modulation r(t) = (1 — cos(27rf/to))/2 is chosen in order 
to minimize the acoustic waves produced during loading. 
The initial void ratio is around v = 0.22. 

The components of the stress are represented by p — 
{<j\ + o 2 )/2 and q ~ (<J\ — 02)/2, where o\ and a 2 are 
the eigenvalues of the averaged stress tensor on the RVE. 
Equivalently, the coordinates of the strain are given by 
the sum 7 — e 2 — £1 and the difference e = —e\ — e 2 of the 
eigenvalues of the strain tensor. We use the convention 
that e > means compression of the sample. The diam- 
eter of the RVE is taken d = 16£ . All the calculations 
were taken in the quasistatic regime. 




FIG. 6. Polygonal assembly confined by a rectangular 
frame of walls. 



A. Superposition principle 

We start exploring the validity of the superposition 
principle presented in Subsec. IIIC. The part (a) of Fig. 



7 shows the loading path during the proportional load- 
ing under constant lateral pressure. This path is also 
decomposed into pieces divided into two parts: one is an 
incremental isotropic loading (Ap = Act and Aq = 0), 
the other is an incremental pure shear loading (Aq = Act 
and Ap = 0). The length of the steps Act varies from 
to OApo to O.OOlpo, where p = 640fcPa. The part (b) 
of Fig. 7 shows that as the steps decrease, the strain 
response converges to the response of the proportional 
loading. In order to verify this convergence, we calculate 
the difference between the strain response of the stair-like 
path 7(e) and the proportional loading path 70(e) as: 

Ae = max |7(e) - 7o(e)|, (46) 




0.2 0.4 0.6 0.8 



e (%) 

FIG. 7. Comparison between numerical responses obtained 
from MD simulations of a rectilinear proportional loading 
(with constant lateral pressure) and stair-like paths. 
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for different steps sizes. This is shown in Fig. 8 for seven 
different polygonal assemblies. The linear interpolation 
of this data intersects the vertical axis at 3 x 10~ 7 . Since 
this value is below the error given by the quasi-static ap- 
proximation, which is 3 x 10~ 4 , the results suggest that 
the responses converge to that one of the proportional 
load. Therefore we find that within our error bars the 
superposition principle is valid. 

A close inspection of the incremental response will 
show that the validity of the superposition principle is 
linked to the existence of tensorial zones in the incre- 
mental stress space. Before this, a short introduction to 
the strain envelope responses follows. 



Ae = 0.0207 Aa /p - 3.1 3x1 0" 7 
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FIG. 8. Distance between the response of the stair-like path 
and the proportional path. 



B. Incremental response 

A graphical illustration of the particular features of 
the constitutive models can be given by employing the 
so-called response envelopes. They were introduced by 
Gudehus [17] as a useful tool to visualize the properties 
of a given incremental constitutive equation. The strain 
envelope response is defined as the image {de = G(da, a)} 
in the strain space of the unit sphere in the stress space, 
which becomes a potato-like surface in the strain space. 

In practice, the determination of the stress envelope 
responses is difficult because it requires one to prepare 
many samples with identical material properties. Numer- 
ical simulations result as an alternative to the solution of 
this problem. They allow one to create clones of the same 
sample, and to perform different loading histories in each 
one of them. 

In the case of a plane strain tests, where there is no 
deformation in one of the spatial directions, the strain en- 
velope response can be represented in a plane. According 
to Eq. (36), this response results in a rotated, translated 



ellipse in the hypoplastic theory, whereas it is given by a 
continuous curve consisting of two pieces of ellipses in the 
elasto-plastic theory, as result from Eq. (35). It is then 
of obvious interest to compare these predictions with the 
stress envelope response of the experimental tests. 

Fig. 9 shows the typical strain response resulting from 
the different stress controlled loading in a dense packing 
of polygons. Each point is obtained from the strain re- 
sponse in a specific direction of the stress space, with the 
same stress amplitude of 10 _4 po- We take qo = 0.45po 
and po — 160kPa In this calculation. The best fit of 
these results with the envelopes response of the elasto- 
plasticity (two pieces of ellipses). For comparison the 
hypoplasticity (one ellipse) is also shown in Fig. 9. 

From these results we conclude that the elasto-plastic 
theory is more accurate in describing the incremental re- 
sponse of our model. One can draw to the same con- 
clusion taking different strain values with different initial 
stress values [38] . These results have shown that the in- 
cremental response can be accurately described using the 
elasto-plastic relation of Eq. (35). The validity of this 
equation is supported by the existence of a well defined 
flow rule for each stress state. 



,x 10" 




-10 -5 5 

e xio" 7 

FIG. 9. Numerical calculation of the incremental strain re- 
sponse. The dots are the numerical results. The solid curve 
represents the fit to the elasto-plastic theory. The dashed 
curve is the hypoplastic fit. 



C. Yield function 

In Subsec. IIIB, we showed that the yield surface is 
an essential element in the formulation of the Drucker- 
Prager theory. This surface encloses a hypothetical re- 
gion in the stress space where only elastic deformations 
are possible [18]. The determination of such a yield sur- 
face is essential to determine the dependence of the strain 
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response on the history of the deformation. 

We attempt to detect the yield surface by using a stan- 
dard procedure proposed in experiments with sand [39]. 
Fig. 10 shows this procedure. Initially the sample is sub- 
jected to an isotropic pressure. Then the sample is loaded 
in the axial direction until it reaches the yield-stress state 
with pressure p and deviatoric stress q. Since plastic de- 
formation is found at this stress value, the point (p, q) 
can be considered as a classical yield point. Then, the 
Drucker-Prager theory assumes the existence of a yield 
surface containing this point. In order to explore the 
yield surface, the sample is unloaded in the axial direc- 
tion until it reaches the stress point with pressure p — Ap 
and deviatoric stress q — Ap inside the elastic regime. 
Then the yield surface is constructed by re-loading in 
different directions in the stress space. In each direction, 
the new yield point must be detected by a sharp change 
of the slope in the stress-strain curve, indicating plastic 
deformations. 

Fig. 11 shows the strain response taking different load 
directions in the same sample. The initial stress of the 
sample is given by qo = 0.5po and Po = 160fcPa. If the 
direction of the reload path is the same as that of the orig- 
inal load (45°), we observe a sharp decrease of stiffness 
when the load point reaches the initial yield point, which 
corresponds to the origin in Fig. 11. However, if one takes 
a direction of re-loading different from 45°, the decrease 
of the stiffness with the loading becomes smooth. Since 
there is no straightforward way to identify those points 
where the yielding begins, the yield function, as it was 
introduced by Drucker and Prager [18] in order to de- 
scribe a sharp transition between the elastic and plastic 
regions, is not consistent with our results. 

Experimental studies on dry sand seem to show that 
the truly elastic region is probably extremely small [4]. 
Moreover, a micro-mechanic investigation of the mechan- 
ical response of granular ratcheting under cyclic load- 
ing has shown that any load involves sliding contacts, 
and hence, plastic deformation [40]. These studies draw 
to the conclusion that the elastic region, used in the 
Drucker-Prager theory to give a dependence of the re- 
sponse on recent history, is not a necessary feature of 
granular materials. 

A question that naturally rises is if the hypoplastic 
theory is more appropriate than the elasto-plastic mod- 
els to describe soil plasticity. Since these models do not 
introduce any elastic regime, they seem to provide a good 
alternative. However, the modern versions of hypoplas- 
ticity depart from the superposition principle, which is 
not consistent with our results. An alternative approach 
to hypoplasticity can be reached from the bounding sur- 
face theory, by shrinking the elastic regime to the cur- 
rent stress point [41]. With this limit one can reproduce 
the observed continuous transition from the elastic to the 
elasto-plastic behavior and in the same time keep the ten- 
sorial zones. However, it has been shown that this limit 



leads to a constitutive relation in terms of some internal 
variables, which lack of physical meaning in this theory. 
In the author's opinion, the necessity to provide a micro- 
mechanical interpretation of these internal variables will 
be important to capture this essential feature of mechan- 
ics of granular materials, that any loading involves plastic 
deformation. 
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FIG. 10. Method to obtain the yield 
surface. Load-unload-reload tests are performed taking dif- 
ferent directions in the reload path. In each direction, the 
point of the reload path where the yielding begins is marked. 
The yield function is constructed by connecting these points. 
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FIG. 11. Strain responses according to Fig. 10. The solid 
lines show the strain response from different reload directions. 
The dashed contours represent the strain envelope responses 
for different stress increments |A<r|. 



VI. INSTABILITIES 

Instability has been one of the classical topics of soil 
mechanics. Localization from a previously homogeneous 
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deformation to a narrow zone of intense shear is a com- 
mon mode of failure of soils [20,42,37]. The Mohr- 
Coulomb criterion is typically used to understand the 
principal features of the localization. This criterion was 
improved by the Drucker condition, based on the hypoth- 
esis of the normality, which results in a plastic flow per- 
pendicular to the yield surface [18]. This theory predicts 
that the instability appears when the stress of the sample 
reaches the plastic limit surface. This surface is given by 
the stress states where the plastic deformation becomes 
infinite. In previous work, it is shown that the normality 
postulate is not fulfilled in the incremental response of 
this model, because the flow and yield directions of Eq. 
(34) are different [38]. Thus, it is interesting to see if the 
Drucker stability criterion is still valid. 

According to the flow rule of Eq. (34), the plastic 
limit surface can be found by looking for the stress val- 
ues where the plastic modulus vanishes. The dependence 
of this modulus on the stress fits to the following power 
law relation [38]: 



Vde, da ■ de > 0. 



(50) 



h = h 



1 



<?o P 



(47) 



This is given in terms of the four parameters: The plas- 
tic modulus ho = 14.5 ± 0.05 at q = 0, the constant 
qo = 0.85 ± 0.05, and the exponents r\ = 2.7 ± 0.04 and 
$ = 0.99 ± 0.02. Then, the plastic limit surface is given 
by the stress states with zero plastic modulus: 

i? 

(48) 



On the other hand, the failure surface can be defined 
by the limit of the stress values where the material be- 
comes unstable. It has been shown that this is given by 
the following relation [38] 



(49) 



q_ 



Here po — l.OMPa is the reference pressure, and q c = 
0.78±0.03MPa. The power law dependence on the pres- 
sure, with exponent (3 = 0.92 ± 0.02 implies a small devi- 
ation from the Mohr-Coulomb theory where the relation 
is strictly linear. 

By comparing Eq. (49) to Eq. (48) one observes that 
during loading the instabilities appear before reaching 
the plastic limit surface. Theoretical studies have also 
shown that in the case of non-associated materials, (i.e. 
where flow direction does not agree with the yield di- 
rection) the instabilities can appear strictly inside of the 
plastic limit surface [16]. In this context, the question 
of instability must be reconsidered beyond the Drucker 
condition. 

The stability for non-associated elasto-plastic materi- 
als has been investigated by Hill, who established the 
following sufficient stability criterion [43]. 



The analysis of this criterion of stability will be pre- 
sented here based on the constitutive relation given by 
Eq. (35). The stability condition of this constitutive re- 
lation can be evaluated by introducing the normalized 
second order work [16]: 



d 2 W = 



da ■ de 



(51) 



Then, the Hill condition of stability can be obtained 
by replacing Eq. (35) in this expression. This results in 

d 2 W = aD-'a + (C ° S( ^ + 0)) cos(fl + > 0, (52) 

where a was defined in Eq. (24). In the case where the 
Drucker normality postulate ip — <fr is valid, Eq. (52) 
is strictly positive and, therefore, this stability condition 
would be valid for all the stress states inside the plastic 
limit surface . On the contrary, for a non-associated flow 
rule as in our model, the second order work is not strictly 
positive for all the load directions, and it can take zero 
values inside the plastic limit surface (i.e. during the 
hardening regime where h > 0). 

To analyze this instability, we adopt a circular repre- 
sentation of d 2 W shown in Fig. 12. The dashed circles 
in these figures represent those regions where d 2 W < 0. 
For stress ratios below q/p = 0.7 we found that the sec- 
ond order work is strictly positive. Thus, according to 
the Hill stability condition, this region corresponds to 
stable states. For the stress ratio q/p = 0.8, the second 
order work becomes negative between 27° < 9 < 36° and 
206° < 6 < 225°. This leads to a domain of the stress 
space strictly inside the plastic limit surface where the 
Hill condition of stability is not fulfilled, and therefore 
the material is potentially unstable. 

Numerical simulations of biaxial tests show that strain 
localization is the most typical mode of failure [7,44]. 
The fact that it appears before the sample reaches the 
plastic limit surface suggests that the appearance of this 
instability is not completely determined by the current 
macroscopic stress of the material, as it is predicted by 
the Drucker-Prager theory. More recent analytic [45] and 
experimental [37,36] works have focused on the role of the 
micro-structure on the localized instabilities. Frictional 
slips at the particles have been used to define additional 
degrees of freedom [45]. The introduction of the par- 
ticle diameter in the constitutive relations results in a 
correct prediction of the shear band thickness. The new 
degrees of freedom of these constitutive models are still 
not clearly connected to the micro-mechanical variables 
of the grains, but with the development of numerical sim- 
ulations this aspect can be better understood. 
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FIG. 12. The solid lines show the second order work as a function of the direction of load for three different stress ratios 
q/p = 0.6 (left), 0.7 (center), and 0.65 (right) with pressure p — IQOKPa. The dashed circles enclose the region where d 2 W < 0. 



VII. CONCLUDING REMARKS 

In this paper we have obtained explicit expressions for 
the averaged stress and strain tensors over a RVE, in 
terms of the micro-mechanical variables, contact forces 
and the individual displacements and rotations of the 
grains. 

A short review on the incremental non-linear models 
has been presented. We emphasize the existence of the 
elastic regime, and the two tensorial zones as predicts the 
theory of elasto-plasticity. We consider also the superpo- 
sition principle of soil mechanics, which is not satisfied 
in the incremental non-linear models. These assumptions 
have been studied using molecular dynamics simulations 
on a polygonal packing. The results are summarized as 
follows: 

• The elasto-plastic theory, with two tensorial zones, 
provided a more accurate description of the incre- 
mental response than the hypoplastic theory. 

• In contradiction to the incremental non-linear mod- 
els, the simulation results show that the superposi- 
tion principle is accurately satisfied. 

• The experimental method proposed by Tatsouka 
has been implemented to identify the yield surface. 
The resulting strain response shows that the tran- 
sition from elasticity to elasto-plasticity is not as 
sharp as the Drucker-Prager theory predicts, but a 
smooth transition occurs. The fact that there is no 
purely elastic regime leads to the open question of 
how to determine the dependence of the response 
of soils on the history of the deformation. 

• The calculation of the plastic limit condition and 
the failure surface shows that the failure appears 
during the hardening regime h > 0. This result is 
analyzed using the Hill condition of stability, which 
states that for non-associated materials the insta- 
bilities can appear before the plastic limit surface. 



These conclusions appear to contradict both the 
Drucker-Prager theory and the hypoplastic models. In 
future work, it would be important to revisit the ques- 
tion of the incremental non-linearity of soils from micro- 
mechanical considerations. 



ACKNOWLEDGMENTS 



We thank F. Darve, Y. Kishino, D. Kolymbas, F. 
Calvetti, Y.F. Dafalias S. McNamara and R. Cham- 
bon for helpful discussions and acknowledge the support 
of the Deutsche Forschungsgemeinschaft within the re- 
search group Modellierung kohdsiver Reibungsmaterialen 
and the European DIGA project HPRN-CT-2002-00220. 



[1] L.D. Landau and E. M. Lifshitz. Theory of Elasticity. 
Pergamon Press, Moscou, 1986. Volume 7 of Course of 
Theoretical Physics. 

[2] P. A. Vermeer. Non-associated plasticity for soils, con- 
crete and rock. In Physics of dry granular media - NATO 
ASI Series E350, page 163, Dordrecht, 1998. Kluwer Aca- 
demic Publishers. 

[3] K. H. Roscoe and J. B. Burland. On the generalized 
stress-strain behavior of 'wet' clay. In Engineering Plas- 
ticity, pages 535-609, Cambridge, 1968. Cambridge Uni- 
versity Press. 

[4] G. Gudehus, F. Darve, and I. Vardoulakis. Constitutive 
Relations of soils. Balkema, Rotterdam, 1984. 

[5] P. A. Cundall and O. D. L. Strack. A discrete numerical 
model for granular assemblages. Geotechnique, 29:47-65, 
1979. 

[6] K. Bagi. Stress and strain in granular assemblies. Mech. 
of Materials, 22:165-177, 1996. 



14 



[7] P. A. Cundall, A. Drescher, and O. D. L. Strack. Numer- 
ical experiments on granular assemblies; measurements 
and observations. In IUTAM Conference on Deformation 
and Failure of Granular Materials, pages 355-370, Delft, 
1982. Balkcma, Rotterdam. 

[8] C. Goldenberg and I. Goldhirsch. Force chains, mi- 
croclasticity, and macroelasticity. Phys. Rev. Lett., 
89(8):084302, 2002. 

[9] K. Bagi. Microstructural stress tensor of granular assem- 
blies with volume forces. J. Appl. Mech., 66:934-936, 
1999. 

[10] M. Latzel. From discontinuous models towards a contin- 
uum description of granular media. PhD thesis, Univer- 
sity Stuttgart, 2002. 

[11] J. P. Bardet. Numerical simulations of the incremental 
responses of idealized granular materials. Int. J. Plastic- 
ity, 10:879-908, 1994. 

[12] Y. Kishino. On the incremental nonlinearity observed in 
a numerical model for granular media. Italian Geotech- 
nical Journal, 3:3-12, 203. 

[13] F. Calvetti, G. Viggiani, and C. Tamagnini. Microme- 
chanical inspection of constitutive modelling. In Consti- 
tutive modelling and analysis of boundary value problems 
in Geotechnical Engineering, pages 187-216., Benevento, 
2003. Hevelius Edizioni. 

[14] M. Oda and K. Iwashita. Study on couple stress and shear 
band development in granular media based on numer- 
ical simulation analyses. Int. J. of Enginering Science, 
38:1713-1740, 2000. 

[15] R. J. Bathurst and L. Rothenburg. Micromechanical as- 
pects of isotropic granular assemblies with linear contact 
interactions. J. Appl. Mech., 55:17-23, 1988. 

[16] F. Darve and F. Laouafa. Instabilities in granular materi- 
als and application to landslides. Mechanics of Cohesive- 
Frictional Materials, 5:627-652, 2000. 

[17] G. Gudehus. A comparison of some constitutive laws for 
soils under radially symmetric loading and unloading. 
Can. Geotech. J., 20:502-516, 1979. 

[18] D.C. Drucker and W. Prager. Soil mechanics and plastic 
analysis of limit design. Q. Appl. Math., 10(2):157-165, 
1952. 

[19] R. Nova and D. Wood. A constitutive model for sand in 
triaxial compression. Int. J. Num. Anal. Meth. Geomech., 
3:277-299, 1979. 

[20] K. H. Roscoe. The influence of the strains in soil mechan- 
ics. Geotechmque, 20(2): 129-170, 1970. 

[21] H. B. Poorooshasb, I. Holubec, and A. N. Sherbourne. 
Yielding and flow of sand in triaxial compression. Can. 
Geotech. J., 4(4):277-398, 1967. 

[22] E. M. Wood. Soil Mechanics-transient and cyclic loads. 
Chichester, 1982. 

[23] Y. F. Dafalias and E. P. Popov. A model of non-linearly 
hardening material for complex loading. Acta Mechanica, 
21:173-192, 1975. 

[24] D. Kolymbas. An outline of hypoplasticity. Arch. Appl. 
Mech., 61:143-151, 1991. 

[25] F. Darve, E. Flavigny, and M. Meghachou. Yield surfaces 
and principle of superposition: revisit through incre- 
mentally non-linear constitutive relations. International 
Journal of Plasticity, 11(8):927, 1995. 



[26] R. Chambon, J. Desrues, W. Hammad, and R. Charlier. 
CLoE, a new rate type constitutive model for geomateri- 
als. Theoretical basis and implementation. Int. J. Anal. 
Meth. Geomech., 18:253-278, 1994. 

[27] W. Wu, E. Bauer, and D. Kolymbas. Hypoplastic con- 
stitutive model with critical state for granular materials. 
Mech. Matter., 23:45-69, 1996. 

[28] I. Herle and G. Gudehus. Determination of parameters of 
a hypoplastic constitutive model from properties of grain 
assemblies. Mechanics of Cohesive-Frictional Materials, 
4:461-486, 1999. 

[29] D. Kolymbas. Modern Approaches to Plasticity. Elsevier, 
1993. 

[30] F. Kun and H. J. Herrmann. Transition from damage 
to fragmentation in collision of solids. Phys. Rev. E, 
59(3):2623-2632, 1999. 

[31] C. Moukarzel and H. J. Herrmann. A vectorizable ran- 
dom lattice. Journal of Statistical Physics, 68(5/6):911- 
923, 1992. 

[32] A. Okabe, B. Boots, and K. Sugihara. Spatial Tessella- 
tions. Concepts and Applications of Voronoi Diagrams. 
John Wiley & Sons, Chichester, 1992. Wiley Series in 
probability and Mathematical Statistics. 

[33] H. J. Tillemans and H. J. Herrmann. Simulating deforma- 
tions of granular solids under shear. Physica A, 217:261- 
288, 1995. 

[34] M. P. Allen and D. J. Tildesley. Computer Simulation of 
Liquids. Oxford University Press, Oxford, 1987. 

[35] E. Buckingham. On physically similar systems: Illustra- 
tions of the use of dimensional equations. Phys. Rev., 
4:345-376, 1914. 

[36] T. Marcher and P. A. Vermeer. Macromodelling of soften- 
ing in non-cohesive soils. In Continuous and Discontin- 
uous Modelling of Cohesive Frictional Materials, pages 
89-110, Berlin, 2001. Springer. 

[37] J. Desrues. Localisation de la deformation plastique dans 
les materieux granulaires. PhD thesis, University of 
Grenoble, 1984. 

[38] F. Alonso-Marroquin and H.J. Herrmann. Calculation of 
the incremental stress-strain relation of a polygonal pack- 
ing. Phys. Rev. E, 66:021301, 2002. cond-mat/0203476. 

[39] Tatsouka F and K. Ishihara. Yielding of sand in triaxial 
compression. Soils and Fundations, 14(2):63-76, 1974. 

[40] F. Alonso-Marroquin and H.J. Herrmann. Ratcheting of 
granular materials. Phys. Rev. Lett., 92(5):054301, 2004. 

[41] Y. F. Dafalias. Bounding surface plasticity. I: Mathemat- 
ical foundation and hypoplassticity. J. of Engng. Mech, 
112(9):966-987, 1986. 

[42] P. A. Vermeer. A five-constant model unifying well- 
established concepts. In Constitutive Relations of soils, 
pages 175-197, Rotterdam, 1984. Balkema. 

[43] R. Hill. A general theory of uniqueness and stability in 
elastic-plas tic solids. Journal of Geotechnical Engineer- 
ing, 6:239-249, 1958. 

[44] J. A. Astr0m, H.J.Herrmann, and J. Timonen. Granular 
packings and fault zones. Phys. Rev. Lett., 84:4638-4641, 
2000. 

[45] H.-B. Muhlhaus and I. Vardoulakis. The thickness 
of shear bands in granular materials. Geotechmque, 
(37):271-283, 1987. 



15 



